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Abstract 

More than half of the sources identified by recent ra- 
dio sky surveys have not been detected by wide-field 
optical surveys. We present a study based on our co- 
added image stacking technique, in which our aim is to 
detect the optical emission from unresolved, isolated 
radio sources of the Very Large Array (VLA) Faint 
Images of the Radio Sky at Twenty-cm (FIRST) sur- 
vey that have no identified optical counterparts in the 
Sloan Digital Sky Survey (SDSS) Stripe 82 co-added 
data set. From the FIRST catalogue, 2116 such radio 
point sources were selected, and cut-out images, cen- 
tred on the FIRST coordinates, were generated from 
the Stripe 82 images. The already co-added cut-outs 
were stacked once again to obtain images of high signal- 
to-noise ratio, in the hope that optical emission from 
the radio sources would become detectable. Multi- 
ple stacks were generated, based on the radio lumi- 
nosity of the point sources. The resulting stacked im- 
ages show central peaks similar to point sources. The 
peaks have very red colours with steep optical spectral 
energy distributions. We have found that the optical 
spectral index a u falls in the range —2.9 < a v < —2.2 
(S v oc v a "), depending only weakly on the radio flux. 
The total integration times of the stacks are between 
270 and 300 h, and the corresponding 5er detection 
limit is estimated to be about m r ~ 26.6 mag. We 
argue that the detected light is mainly from the cen- 
tral regions of dust-reddened Type 1 active galactic nu- 
clei. Dust-reddened quasars might represent an early 
phase of quasar evolution, and thus they can also give 
us an insight into the formation of massive galaxies. 
The data used in the paper are available on-line at 
http : //www. vo . elte .hu/doublestacking. 

1 Introduction 

A significant fraction of sources that have been de- 
tected in recent deep radio sky surveys have not yet 
been detected in the optical wavelength regime by wide 
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field surveys. From the about 10 6 objects identified by 
the Very Large Array (VLA) Faint Ima ges of the Ra- 
dio Sky at Twenty - cm (F IRST) survey 1 Becker et all 



1995 IWhite et all Il997ft . only approximately 30 per 



cent have optical count erparts in the legacy Sloan Dig- 
ital Sky Survey (SDSS; lYork et all l2000l: llvezic et al 
l2002h at a limiting magnitude of r ~ 22.2. We estimate 
that for the SDSS Stri pe 82 co-added survey , which is 
about 1.8 mag deeper ( Abazaiian et al. . 20091 ). this de- 
tection ratio is about 42 per cent. 

The nature of the optically undetected frac- 
tion has been the subject of debate for many 
years. Recent studies have shown that there ex- 
ists a sign ificant populati o n of re d dened radio- se lected 



quasars (IWebster et all Il995l: ICutri et all 12001 



2003 



2005 



Glikman et al 
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2004; 
20071) . 
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ally accepted unification model (e.g. 
a likely explanation for the low optical luminosity of 
these quasars is that their central engines are obscured 
by an optically thick dust torus. 

In this paper, we use a new image stacking tech- 
nique to investigate the average optical properties of 
unresolved faint point sources from the VLA FIRST 
survey that have no detected optical counterparts in 
the SPSS Stripe 82 (equat orial stripe) co-added data 
set ( Abazaiian et al. , 2009) We call our technique deep 
co-add stacking (DCS) because we stack images that 
have already been co-added from multiple observations 
of the SDSS equatorial stripe. In our technique, we 
put strong emphasis on accurate background estima- 
tion and on correcting for selection effects. 



1.1 Image stacking 

Co-adding repeated observations and stacking of im- 
ages of different, but similar objects can be used to 
extend the limiting magnitude and surface brightness 
limits of surveys. When stacking images of different ex- 
tended sources with varying apparent sizes - depend- 
ing on the specific needs - cut-out images are made, 
which can be rotated, resized and/or otherwise scaled 
together to make the individual objects overlap. Stack- 
ing increases the signal-to-noise ratio of the images by 
reducing the photon noise. For undetected sources, 
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stacking can be successfully used to reveal an average 
image of extremely faint sources that otherwise would 
be beyond the detection limits of the instruments. The 
individual objects must be similar to each other in or- 
der to be able to give a physical interpretation to the 
average images. The application of image stacking is 
limited by a couple of factors, of which the thermal 
and read-out noise of the detectors are the most impor- 
tant. However, modern detector technology makes it 
possible to stack hundreds of images together, without 
significant systematic background noise in the results. 

One important disadvantage of image stacking is its 
limited capability for statistical analysis. In the op- 
tical wavelength regime, usually only one or a few 
stacked images per band are obtained from the indi- 
vidual cut-outs. Bootstrapping techniques can be used 
to determine the variance of the measured properties of 
the stacked images by sacrificing the surface brightness 
depth and signal-to-noise ratio. 

When we are interested in faint objects only, it is 
important to exclude those pixels of the individual cut- 
outs from stacking that are brighter than a carefully de- 
termined threshold; otherwise, the strong signal from 
the brighter sources or stray light from nearby bright 
sources would wash away the very weak signal of the 
faint object. Fine-tuned masking techniques are re- 
quired to address this problem, as we show in Sec- 
tion E2 

Previously, image stacking has been used success- 
fully to recover the undetected light from objects that 
are visible in a given wavelength range but are too 
faint i n another to detect them directly. IWhite etal 



(|2007D stacked VLA FIRST images of optically identi- 
fied quasars, and investigated the radio properties (the 
correlation between radio an d optical lum i nositi es, ra- 
dio loudness) of the sample. Hodge et al. ( 20081 ) iden- 
tified the faint radio emission from a diverse sample 
of quiescent (optically unclassifiable) galaxies from the 
SDSS by stacking their FIRST radio images. They 
showed that the radio emission from these galaxies can 
be related to their star formation rate, or a signifi- 
cant part of these galaxies might harbour quiescent 
active galactic nuclei (AGNs) in their cores and their 
activity is related to the stellar mass. They continued 
their analysis by stacking radio images of luminous red 
galaxies (LRGs) from the SDSS. They showed that low 
radio luminosity AGNs (with 1400 MHz flux densities 
in the 10 < >Sint,i400 < 100 juJy range) are characteris- 
tic of the LRG population, and that an evolution of the 
nuclear activity is apparent for red shifts in the range 
0.45 < z < 0.6 (iHodee et all l2009h . Radio contribu- 



tion by star formation was ruled out by sample selec- 
tion (i.e. LRGs are passive galaxies with no significant 
star formation). 

Another use of the stacking met hod is the detection 
of faint haloes of extended objects. IZibetti et all (|2004h 
stacked SDSS optical images of edge-on disc galax- 
ies to detect their stellar halo components, and they 



reached a surface brightness level as low as jjL r ~ 31 
mag arcsec -2 . Later, they analysed the radial pro- 
file of intracluster light by stacking the SDSS images 
of 683 galaxy clusters. They identified faint emission 
from intergalactic stars as far as 700 kpc from the 
centres of the clusters, by pushing the surface bright- 
ness limit of the stacked images to n r ~ 27.5 — 30 



mag arcsec 2 in the SDSS r-band (IZibetti et alll2005l) 



Also, IZibetti et al.1 (]2007l ) applied their technique to 



reveal light from Mgn absorbers. lHathi et al.1 (^2008) 
used image stacking to determine the average prop- 
erties of the surface brightness profiles of very dis- 
tant, compact gala xies in the Hubble Deep Field. 
Bergvall et al. I (l2010h detected the extended red haloes 
of low surface brigh t ness S PSS galaxies. Most recently, 
Tal Sz van Dokkum d201lh have investigated the faint 
stellar haloes of LRGs by stacking SDSS images. They 
have found that stellar light in massive elliptical galax- 
ies can be traced out to about a radius of 100 kpc. 

An interesting similarity of these SDSS-based studies 
is that the st acked images show e xtremely red colours. 
According to Zibetti et al.l ( 20041 ) the r — i colours of 
the stellar haloes of the galaxies reach up to about 
0.8 mag, which is attributed to very old and/or metal- 
abundant stellar populations. As an explanation , for 
the case of edge-on disc galaxies, |de Jond (|2008h has 
provided evidence that the anomalous halo colours are 
at least partly a result of the underestimated effect of 
the extended tails of the point spread function (PSF). 



1.2 Radio— optical cross-identification 

Cross-identification of radio sources with their optical 
counterparts is a challenging problem because of the 
complex morphology of radio galaxies. Many studies 
have previously investigated the possibilities of auto- 
matic radio-optical cross-identification and the connec- 
tions between th e optical and radio prop erties of extra- 
galactic objects. iMcMahon et al. ( 20021 ) automatically 
identified the optical counterparts of the FIRST radio 
sources in the Cambridge Automated Plate Measure- 
ment (APM) scans of the First Palomar Observatory 
Sky Survey (POSS-I) plates, and they found about 
70.000 matches. According to them, in the case of 
the aforementioned catalogues, positional coincidence 
is enough to cross-identify radio point sources with op- 
tical sources at false-positive levels lower than 5 per 
cent when an association radius of 2 arcsec is used. 
This false-positive rate is, of course, only for the par- 
ticular combination of catalogues of that study. The 
rate depends on the astrometrical acc uracy of the cat 



alogue s and the density of sources. IMcMahon et al 
(|2002l ) also addressed the cross-identification of double 
radio sources with their optically observable counter- 
parts. They found that an optical counterpart is likely 
to be found halfway between the components of the 
radio source pairs. 
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1.3 Origins of the unresolved radio 
emission 

In Section [21 we construct a sample of optically un- 
detected, unresolved radio sources with integrated 
1400 MHz radio flux above 6*1400 > 1 m Jy- We con- 
sider several types of radio sources that can appear un- 
resolved in the VLA FIRST images: galactic pulsars, 
radio stars, hotspots of radio jets associated with the 
outer regions of radio galaxies and AGNs (including 
quasars). 

A flux density limit of 1 mJy is useful to separate 
AGNs from low redshift starburst galaxies. This limit 
is sufficiently bright to exclude most of the starburst 
galaxies, and yet faint enough to pe rmit high redshift 



AGN s and classical r adio g alaxies ([Windhorst et al 
19851 : iHopkins etail 12000) Th e 1-mJy limit was 



used by Waddington et al.T ( 200ll) . who compiled the 
Leiden-Berkeley Deep Survey (LBDS) Hercules sam- 
ple, which is still the largest optically complete sample 
of mJy radio sources at 1400 MHz. 

Because we restrict our analysis to the area covered 
by the SDSS Stripe 82 (-50° < a < 59°, \S\ < 1.26°), 
the high galactic latitude and the relatively high se- 
lection limit on the radio flux make it very unlikely 
that the resulting point sources are pulsars. For exam- 
ple, the comprehensive Australia Telescope National 
Facility (ATNF) pulsar catalogue ( Manchester et aD . 
120051) contains only two pulsars within the footprint of 
Stripe 82. This number is negligible compared to the 
more than 20.000 unresolved radio objects that are cat- 
alogued in the VLA FIRST data set within the same 

footprint. 

~ (fl999l) . 



According to iHelfand et al 



there are 26 

known radio stars brighter than 0.7 mJy at 1400 MHz 
discovered in the FIRST survey, from which only three 
sources lie in the Stripe 82 footprint. Therefore, it is 
very unlikely that our optically undetected sample of 
2116 objects contains any radio stars. 

Hotspots of radio galaxy jets might also appear as 
point sources in FIRST images. They are almost al- 
ways associated with an easily detectable nearby ellip- 
tical galaxy and they tend to appear in pairs, on both 
sides of the galaxies. We have identified many sources 
that are very good candidates for such hotspots, but 
we have excluded them from the present analysis using 
a method described in Section [2TTTI 

The fourth type of objects that appear as unresolved, 
isolated radio sources are the central regions of super- 
massive black holes currently accreting matter. In Sec- 
tion [HI we argue that we detect the optical light from 
this type of sources. 

1.4 Reddened radio-loud quasars 

According to the widely accept ed AGN u nifica- 
tion scheme (fo r a r eview, see lAntonuccll . Il993 : 



properties. Type 1 AGNs have both broad and narrow 
emission lines and show blue optical continua. This 
implies that the central regions of these objects are ob- 
served directly with no significant quantity of obscur- 
ing material along the light of sight. However, Type 2 
sources show red optical spectra and narrow emission 
lines only. It is believed that Type 2 AGNs are viewed 
from such directions that the dense dust clouds sur- 
rounding the central engines fall into the line of sight 
and obscure much of the visible light coming from the 
accretion discs. The broad emission lines are associated 
with the central regions of the accretion discs, where 
high velocity motion is responsible for the broadening 
of the lines. Polarimetric studies have revealed that 
weak, highly polarized broad emission lines are also de- 
tectable in the spectra of Type 2 AGNs. These highly 
polarized lines are attributed to the light that origi- 
nates from the central regions but is reflected from gas 
clouds at greater distances from the centres. 

It has been demonstrated previously that there ex- 
ists a strongly reddened population of quasars, which 
are often missed by the conse rvative quasar sele ction 
criteria of optical sky surveys. Ilvezic et al. (2002) pre- 
sented a very comprehensive study about the optical 
and radio properties of positionally matched FIRST 
and SDSS sources. By analysing the distribution of 
quasars, they pointed out that the existence of a sig- 
nificant population of highly obscured quasars that are 
visible in FIRS T but not in SPSS c annot be ruled 
out by the data. iRichards et al. I (120031) estimated that 
about 10 per cent of the dust-reddened quasars are 
missing from the SDSS catalogue because of selec- 
tion effects. They also estimated that the fraction of 
broad absorption line (BAL) quasars rises with red- 
der colour and can reach 20 per cent for reddened 
quasars. In the case of BAL quasars, the BALs are 
most likely to originate from the fast moving dust ex- 
pelled from the central regions by the strong wind 
of th e nucleus (IHazard et all. Il984 IWevmann et al 



199ll ISpravberrv fe Folt zL Il992t lOele et all Il99£ : 



Schmidt fe Hinesl.ll999HBecker et all l2000HHall et al 
20021 ; iTrump et all 1200 



Glikman et all ([20071 ) have confirmed that the selec- 



tion of quasars by infrared colours yields a significantly 
redder quasar population than selection based on op- 
tical surveys. They have also reported, based on two- 
frequency radio observations, that the red colours of 
these quasars are more likely to be a result of extinc- 
tion by dust rather than enhanced synchrotron emis- 
sion. They have found that the intrinsic extinction 
values of these reddened quasars can reach as high as 
E(B - V) ~ 2.5 mag. 

Certain models suggest that AGNs are turned on by 
galaxy mergers and that they spend a significant time 
obscured by the dust that originates from the merg- 



Urrv fe PadovaniL I1995L and references therein) , AGNs 



ing h ost galaxies (j Sanders et all . [1988; Hop kins et al 
2005). In this phase, the AGNs have high intrin 



are divided into two main classes based on their optical sic luminosities, but they appear as heavily reddened 
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Type Is (or ultraluminous infrared galaxies) because of 
strong obscuration by the dust of the hosts. The unob- 
scured quasars become visible once the radiation from 
t he central engines wip es out the obscuring material. 
Urrutia et al. (2009) have investigated the proper- 



ties of heavily reddened quasars of a sample com- 
piled fro m FIRST, the Two- M icron All-Sky Survey 
(2MASS; ISkrutskie et "all . l200fih and SDSS. The ma- 
jority of their sample was spectroscopically confirmed, 
and they found that a quasar population with high in- 
trinsic extinction does exist far from t he traditional 
photo metric quasar selection criteria. lUrrutia et al 



(2009) measured the reddening values to be in the 
range 0.1 < E(B — V) < 1 mag, where the higher val- 
ues are more consistent with obscuration by the host 
galaxies rather than being a result of the dust torus 
surrounding the central engines. They have also found 
that these quasars are likely to show BALs of elements 
at low ionization states. This observation suggests that 
these AGNs are at an early period of their active cy- 
cles, when strong winds driving the dust out from the 
central regions cause the BALs. 

For a comprehensive review of the spectral en- 
ergy distributions (S EDs) of Type 1 quasars, see 



Richards et ail (|2006l ). 



1.5 Structure of the paper 

In Section^ we describe the radio object sample se- 
lection method, and the source of the optical imaging 
data. We give details of the image stacking algorithm 
in Section [31 and we present the results of the analysis 
- including photometry, SED and radial surface bright- 
ness profiles of the stacked objects - in Section [5j Fi- 
nally, in Section [6l we discuss the optical and infrared 
properties of our sample. 

We use the concordance flat A cold dark matter 
(CDM) cosmology throughout this paper, with param- 
eters of H = 71 km s _1 Mpc" 1 , Q M = 0.27 and 
n A = 0.73. 

We have developed our own software tools, written 
in IDL, to generate cut-out images and masks, and for 
the entire image stacking process. The code is available 
from the authors on request. 



2 Data 

For the present analysis, we u se 1400 MHz radio data 
from the VLA FIRST survey (|Becker et all Il995h and 
u, g, r, i and z-band optical images from th e SPS S 
Stripe 82 co-added data set ( Abazaiian et all [2009). 
The aim of our analysis is to detect the faint optical 
emission of sources that are individually undetected in 
the optical, but which have significant emission in the 
radio. We choose our sample selection criteria accord- 
ing to the requirements detailed in the following two 
subsections. 




r [arcsec] 

Figure 1: Pair-correlation function of 6019 FIRST ra- 
dio point sources brighter than <Si n t.i4oo > 1 m Jy ; plot- 
ted for small separations only. The reason behind the 
strong correlation at very small angles is that multi- 
ple point sources are likely to be associated with the 
same radio galaxy. The vertical dashed line indicates 
our 1.5 arcmin limit on the minimum separation from 
any other radio sources. 



2.1 Sample selection from the FIRST 
catalogue 

We only select those FIRST objects that are within 
the SDSS equatorial stripe footprint (-50° < a < 59°, 
\6\ < 1.26°). Objects are required to be compact ra- 
dio sources, i.e. point-like without any extended ra- 
dio features. The FIRST catalogue assigns a side-lobe 
probability Ps to each object based on oblique decision 
tree classifiers (refer to the FIRST web-sit43 for de- 
tails). We require that the side- lobe probability must 
be less than Ps < 0.1. In addition, the integrated 
1400 MHz fluxes of the objects must be higher than 
<Sint,i4oo = 1 mjy. We visually inspect all automati- 
cally selected sources, as described in Section l2~2l 

In Fig. [T] the angular pair-correlation function of the 
FIRST radio point sources is plotted for small sepa- 
rations. A strong correlation at very small angles is 
clearly visible, but the function becomes flat above 
1.5 arcmin. This behaviour originates from the sim- 
ple fact that radio galaxies have complex morpholo- 
gies with many bright features: the multiple lobes and 
hotspots of radio galaxies are likely to appear as dou- 
ble or multiple associations of unresolved sources in 
the FIRST images that are close to each other. In most 
cases, the optical counterparts of these galaxies are eas- 
ily identifiable in the proximity of the radio sources. 
In order to make sure that only isolated sources are 
included in the sample, we consider the angular pair- 
correlation function of them, and further restrict the 
selection conditions to exclude radio sources with an- 
other close radio companion within 1.5 arcmin. This 



: http : //sundog. stsci . edu/ 
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criterion is though to effectively exclude most point 
sources that are associated with optically resolved ra- 
dio galaxies. 

Because we are looking for compact radio sources 
with no clearly identified optical counterparts, in the 
analysis we only include those radio sources that do 
not have any matches in the SDSS Stripe82 co-added 
catalogue within 3 arcsec. This condition is much less 
restrictive than the typica l astrometric error of S DSS, 



which is about 0.1 arcsec (jStoughton et al.. 2002|), but 



it is strict enough to exclude any false positive matches. 
It also excludes extended radio sources with small di- 
ameters that might show an offset of a few arcsec be- 
tween the optical and radio source positions. This 
value is chosen because practically all true SDSS- 
FIRST matches ha ve smaller separations than 3 arcsec 
(Jlvezic et all |2002|) . 
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Table 1: Sources of the final DCS sample. (The 
entire table is available in machine-readable form at 
http://www.vo.elte.hu/doublestacking. A portion is 
shown here.) 



2.2 Visual inspection of the images 

After applying the former criteria, 2626 sources re- 
main. We identified these radio source s in th e Deep 
VLA Stripe 82 catalogue ((Hodge et all l201ll) which 
has double the resolution and about three times the 
signal-to-noise ratio of the FIRST survey. Although 
the deeper survey does not cover the entire Stripe 82 
region (only about 120 square degrees of the total 280 
square degrees, 43 per cent), we still can use the deeper 
catalogue to verify the efficiency of the selection cri- 
teria we used to filter the FIRST data set. We are 
especially interested whether radio point sources can 
be segregated from extended sources based solely on 
the information available in the lower resolution FIRST 
catalogue. Therefore, we visually inspect the FIRST, 
Deep VLA and SDSS Stripe 82 images of all previously 
selected sources. 

First of all, we exclude all radio objects if the cor- 
responding optical frames are visibly affected by scat- 
tered light from nearby bright stars. The number of 
these objects amounts to about 180. In a few cases, we 
have to exclude radio objects because they overlap with 
the outer regions of very bright, extended galaxies. 

About 10 per cent of the automatically selected 
FIRST objects is rejected during visual inspection be- 
cause they clearly show resolved features or appear 
very faint and different from the average FIRST PSF. 

Next, we compare the FIRST and Deep VLA im- 
ages of those objects which were observed during both 
surveys. While all objects appear pointlike in the 
FIRST images (although some show a slight eccen- 
tricity), about 4 per cent of the FIRST point sources 
are clearly resolved into radio galaxy features in the 
deeper and higher resolution VLA survey. Unfortu- 
nately, there appears to be no way to identify the re- 
solved sources based solely on their FIRST catalogue 
properties. Hence, exclusion of radio objects that be- 
come resolved at higher resolution is only possible in 
one half of the sample. 
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Table 2: Distribution of the DCS sources by subsam- 
ples. 5i„ t) i4oo is the 1400 MHz flux density interval of 
the bin; 0-5 is the variance of the flux density, N is 
the number of the sources of each subsample; and t exp 
the total effective exposure time of the stacked optical 
images. 
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Figure 2: Distribution of the radio flux Si n t,i4oo of the 
DCS sample. The grey areas indicate the three sub- 
samples, light grey on the left: 1 — 2 mJy; dark grey: 
2 — 4 mJy; light grey on the right: > 4 mJy. 
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To select a radio source we required that it should 
be at least 1.5 arcmin away from any other FIRST ra- 
dio sources. When looking at the deeper radio survey 
images, objects of low radio brightness are apparent 
within 1.5 arcmin in a few cases. Also, a very few ob- 
jects detected in the FIRST survey are missing from 
the Deep VLA images. We attribute these false detec- 
tions to glitches in the FIRST observations. 

By looking at the higher- quality Deep VLA images, 
it is estimated that about 9 per cent of the FIRST 
objects would be excluded at the depth of the Deep 
VLA observations. However, we have found that the 
exclusion or inclusion of these objects does not mod- 
ify our final findings significantly; the difference in the 
measured fluxes is less than 5 per cent, which is in the 
range of the estimated error value. Although it can 
only be done for half of the data set, we still decide to 
exclude them from our final stacks. 



2.3 Subsample definitions 

The final, deep co-add stack (DCS) sample consists 
of 2116 objects, listed in Table [TJ The sources are 
grouped into three, roughly equal cardinality subsam- 
ples based on their integrated apparent radio flux den- 
sities 5'int,i4oo- The size of the subsamples make it 
possible to generate stacks with high enough signal- 
to-noise ratio. The definitions of the subsamples are 
based on apparent flux because there are no redshift 
estimates available for the sample. Still, we can ex- 
pect to see some dependence of the derived parameters 
on the apparent radio flux at the end of the analysis 
if the objects lie in a relatively narrow redshift range 
(i.e. a correlation between apparent and absolute radio 
luminosities exists). 

Table [5] summarizes the properties of the three DCS 
subsamples. The distribution of the radio fluxes is plot- 
ted in Fig. H 



2.4 SDSS co-added images 

For the stacking, we use the SDSS equatorial stripe co- 
added imaging data to create cut-outs centred on the 
selected DCS source coordinates. The SDSS Stripe 82 
images were co-added from abo ut 20-40 individual ex- 
posures (jAbazaiian et al. I, [2Q09)| Each exposure took 
54 s and the resulting catalogue has an estimated r- 
band magnitude limit of about r < 24 mag, which is 
about one magnit ude deeper than the rest of SDSS 
(|Gunn et all 1 19981) . 



2 The number of scans of Stripe 82 now exceeds 70, but only 
20-40 scans have been included in the co-addition (mostly data 
collected in 2005). 



3 Image stacking 

3.1 Cut-outs 

The DCS sample consists of radio sources with opti- 
cal emissions below the detection limit of the SDSS 
Stripe 82 co-added catalogue. To reveal the average 
undetected light from these sources, first we generate 
cut-out images from the already co-added Stripe 82 im- 
ages centred on the radio coordinates. Cut-outs of the 
size 200 by 200 pixels are made, which is equivalent 
to 80 by 80 arcsec at the resolution of SDSS. Cut-out 
image examples with radio flux contour overlays are 
presented in Fig.S 

Because the background has already been subtracted 
from the images during the co-addition, the cut-out 
images should have sky levels of 0. Later, we show 
that this is not the case. In Section [231 we introduce 
a technique to re-estimate the zero level of the cutouts 
using an algorithm based on using cut-outs taken at 
random coordinates. 

3.2 Masking 

A significant portion of the pixels of the 200 by 200 
Stripe 82 cut-outs is from bright objects that are out 
of our interest. Also, very frequently, stray light from 
nearby, bright (pii < 16 mag) stars contaminates the 
entire image, but light from interstellar clouds of the 
Milky Way can also be seen in some of the cut-outs. 
In order to reveal the very weak signal from the un- 
detected sources, we have to exclude bright pixels effi- 
ciently. Our masking algorithm is based on a histogram 
cutting technique as described below. 

3.2.1 Bright objects 

The 1.5 arcmin limit on the closest radio source im- 
plicitly assures that no associated radio galaxies will 
appear in the images (c.f. Section 12.11) . Also, the 
condition that the selected radio sources have to be 
further away than 3 arcsec from any detected optical 
sources ensures that there will be no readily visible op- 
tical sources in the centres of the cut-out images. How- 
ever, outside the 3 arcsec radius there can be numerous 
objects that have to be masked out before stacking the 
cut-outs together, otherwise the weak signal from the 
centres would not be detectable. 

3.2.2 The masking algorithm 

First, we have to determine the pixel values above 
which pixels should be rejected. Setting the right value 
for the masking threshold is crucial for the stacking to 
be successful. The values should be chosen in such 
a way that all pixels belonging to readily detectable 
objects are masked, while sky pixels and pixels of un- 
detectable faint objects are kept. 

The histogram of the pixel values of a single Stripe 
82 co-added frame is plotted in panel (a) of Fig. @] as 
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Figure 3: Sample cut-out images from the SDSS equatorial stripe (Stripe 82) co-added survey with 1400 MHz 
radio flux contour overlays from FIRST. The top two images are examples of simple, isolated radio point 
sources, while the lower two ones show radio galaxies with more complex radio features. Only the top left 
object contributes our sample, the others (including the top right one) are excluded by the criterion which 
requires that no other radio objects should be in a 1.5 arcmin radius of a given point source. The images are 
80 arcsec wide on the side. The white circles have 3 arcsec radii. The contour levels are 0.4, 0.65, 1.0, 1.45, 2.0, 
2.65, 3.4, 4.25, 5.2, 6.25, 7.4 mJy beam" 1 . 
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Figure 4: Histograms illustrating the steps of our masking algorithm. Note the logarithmic vertical scale, a) 
Histogram of the pixel values of a single, typical SDSS Stripe 82 co-added frame, b) Histogram of the same 
frame convolved with a circular top-hat kernel with a radius of 1 arcsec. c) Histogram of the frame after applying 
the mask (solid line). The unmasked histogram is indicated with a dotted line. The solid vertical lines are at 
the mean, the dashed lines are at the ±1<7 position of the original distribution. Sec the text for more details. 
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Figure 5: Left: a sample cut-out image from the stack. 
Right: The same image overlaid with the mask (shown 
as black regions). 

an example. Because we use co-added images, the sky 
level is already subtracted at least to first order. This 
can be seen clearly from the figure; the histogram peaks 
nearly at zero and background pixels can have both 
positive and negative values. The negative part of the 
histogram is completely dominated by the noise of the 
sky pixels, while the longer tail towards the bright pix- 
els comes from the objects. 

The masking process is carried out in four steps, as 
follows. 

1. First, we determine the masking threshold. We fit 
the negative side of the distribution of pixel val- 
ues with a Gaussian function, assuming that the 
distribution of the background pixels is symmetric 
around zero. The masking threshold is then set to 
be at the la variance of the fitted Gaussian. We 
indicate the value of the mean and ±1<7 as vertical 
lines in the panels of Fig. [U 

A higher threshold would let more bright pixels 
into the statistics, which would increase the contri- 
bution of light from unwanted objects to the back- 
ground, resulting in lower contrast stacks. How- 
ever, if the value of the threshold was too low, we 
would filter out much of the light from the inter- 
esting faint sources themselves and render them 
undetectable in the stacks. We have found that 
the la threshold yields a good balance between 
contrast and detectability. In Section POl we show 
that changing this value slightly does not affect the 
final results. 

The a values for each SDSS band are determined 
from the pixel distribution of the whole ensemble 
of cut-outs, and the same values are used through- 
out the analysis. 

2. Secondly, we make smoothed images from the orig- 
inal cutouts. This step is important because the 
imaged sources have diffuse edges. Therefore, 
masks derived directly from the cut-out images 
by cutting at the masking threshold would not be 
appropriate. The masks have to be slightly ex- 
panded to cover the smooth edges of the objects. 



We also want to avoid the masks being too frag- 
mented. To achieve these goals, we convolve the 
cut-out images with a circular top-hat kernel with 
a radius of 1 arcsec. The histogram of a convolved 
image is plotted in Fig. |U(b). 

3. Thirdly, we create the masks from the convolved 
images; pixels with values above la are masked. 
Without the initial convolution of the original im- 
ages with the top-hat kernel, masked areas would 
be strongly non-contiguous because the values of 
the masking thresholds are chosen to be very close 
to the background noise level. Because the radius 
of the smoothing kernel limits the minimum ra- 
dius of the objects that are masked out, we al- 
low the kernel to have a radius just slightly larger 
than the typical seeing of SDSS. We have found 
that a kernel radius of less than 1 arcsec would 
make the masks too fragmented. The chosen ker- 
nel makes the masked areas more contiguous and 
makes them look more or less circular with sharp 
edges. A typical mask is presented in Fig. [51 

4. Finally, we generate the masked images by multi- 
plying the original (non-convolved) cut-outs with 
the image masks. The histogram of a masked 
image is plotted in Fig. @|c) with a solid line. 
The dotted line shows the original, unmasked his- 
togram. It is clear from the figure that the original 
symmetric distribution of the background is pre- 
served during the masking, while object pixels are 
effectively masked out. 

3.3 Stacking 

We create the stacked images from the co-added Stripe 
82 cut-outs centred on the DCS source coordinates. 
The cut-outs for each SDSS band are grouped together 
into three stacks based on the radio flux density of the 
DCS objects as described in Section [2] First, bright 
pixels are masked out as explained in Section l3~2l then 
the remaining pixel values are averaged. Also, the stan- 
dard deviation of the pixel values at every pixel coor- 
dinate is determined to characterize the error. 

Stacked images are calculated according to the fol- 
lowing formula: 



2^=1 w ikl 

Here Ski is the pixel value of the stack, pm is the pixel 
value of the i th cut-out image and the subscripts k and 
I index the pixel coordinates (k, I = [1,200]), Wiki is 
the mask, its value is if the pixel is masked out and 
1 if it is not; N is the total number of cut-outs. It has 
been mentioned previously that, although the sky level 
was subtracted from the Stripe 82 co-added images, the 
average background levels of the cut-outs usually show 
a small bias. The correction term Cki is introduced to 
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correct for this bias. We will explain this last term in 
detail in Section I3T41 

In Fig. [5J we plot the stacked images for every radio 
flux bin and SDSS imaging band. A PSF-like central 
object is clearly visible in the centres of the images. 
We attribute these peaks to the optical emission from 
the radio-selected objects. 

3.4 Sky level recalibration 



The SDSS imaging pipeline (jStoughton et all . 120021 ) es- 
timates the local sky level using clipped median in 256 
by 256 pixel boxes, centred on every 128 pixels of the 
2048 by 1489 pixel framed. A sky constructed by inter- 
polation of these local sky estimates is then subtracted 
from the frames, prior to co-adding the frames of dif- 
ferent scans. The resolution of the sky estimation does 
not allow for complete elimination of the light from 
high latitude galactic clouds due to their finer texture. 
Although the imaging pipeline attempts to subtract 
the wings of the PSF of saturated stars, the algorithm 
is not too aggressive. Wings of stars can easily fill the 
80 by 80 arcsec cut-outs entirely. These two effects can 
significantly increase the average background level and 
reduce the contrast of the stacked images, thus need to 
be eliminated. 

The cut-outs we generate show a tiny bias of the 
background level towards negative values. Because we 
are dealing with objects that are extremely faint in 
the optical bands (in the range of 30-300 nJy arcsec -2 
or 25-28 mag arcsec -2 ), this small bias (estimated to 
be about 10-60 nJy arcsec -2 or 27-29 mag arcsec -2 , 
depending on the photometric filter used) would sig- 
nificantly affect the photometry of the stacked images. 
The bias is likely to be the consequence of a remaining 
gradient in the background level in the whole co-added 
Stripe 82 frames, which becomes apparent if only small 
parts of the frames are considered. Also, the algorithm 
applied to co-add the Stripe 82 images has used 2.326ct 
clipping to determine the sky level, while we use a much 
restrictive la masking threshold. As a consequence, 
the sky level of the cut-outs has to be re-estimated. 
To do this, we simply subtract the average value of 
the unmasked pixels from the cut-outs on a per image 
basis. 



3.5 Selection bias of the stack sky levels 

When composing the DCS sample, we imposed se- 
lection criteria based on the separation of the radio 
sources from the neighbouring optical sources. This 
introduces a strange bias to the background of the cut- 
outs, as we will show in this section. 

To investigate the issue, we generate stacked images 
of cut-outs centred on random coordinates. We call 
these special stacks random stacks. First, for every 



DCS source coordinate pair, a random coordinate pair 
is generated within the same Stripe 82 frame where 
the DCS object is. We require these random points to 
satisfy the same selection criteria as the DCS objects, 
i.e. they are at least 3 arcsec away from any opti- 
cally detected Stripe 82 co-added sources and at least 
1.5 arcmin away from any other FIRST radio sources. 
Second, random cut-outs are masked as described in 
Section I5~2l Third, the pixel values are averaged ac- 
cording to the following formula: 



'fei 



]U=i win 



(2) 



3 For more details, see 
http: //www. sdss . org/dr7/algorithms/sky . html 



where s]^ denotes the pixel value at the k and I pixel 
coordinates of the random stack image (the upper R is 
used to identify that this is a random stack); pf~ kl is the 
pixel value of the random cut-outs, where i indexes the 
individual cut-outs; Wiki is the mask, its value is is 
the pixel is masked out and 1 if it is not; N is the total 
number of cut-outs, as before. Separate random stacks 
are generated for all five SDSS bands and for all three 
DCS subsamplcs. (Note, that it is important that the 
random coordinates are chosen to be within the same 
Stripe 82 frame as the original DCS source coordinates 
are, so different random stacks are necessary for each 
subsample.) 

To characterize the variance of the random stacks, 
we generate 50 of them using different sets of ran- 
dom coordinates. To get the correction term Cki in- 
troduced in Section 13.31 we combine the 50 random 
stacks into a super-stack. The super-stack (denoted as 
Cm in Eq. [1} is simply the pixel- wise average of the 50 
random stacks. 

The pixels of these final random super-stacks are 
once more averaged in annuli centred on the middle 
of the images to determine the radial surface bright- 
ness profile of the background. The resulting profiles 
are plotted in Fig. [3 

The significant feature of the radial surface bright- 
ness profiles are the dark spots at the centres of the 
images. They are very well visible even though vari- 
ance is higher at these small radii due to the lower pixel 
counts. We account these dark spots to the selection 
criterion that restricts coordinate selection to coordi- 
nates that are at least 3 arcsec away from the centres 
of the closest Stripe 82 objects. The plateaus of the 
surface brightness profiles of the random stacks come 
from the faint light of nearby objects. Contribution 
to these plateaus from stray light from farther bright 
objects is also expected. 

4 The reliability of the stacking 
method 

Simple averaging is known to be sensitive to outliers, 
thus reliable masking of the bright pixels of cut-outs 
prior to averaging is inevitable. We apply three tests 
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Figure 6: The stacked images of the radio-selected cut-out images, for each radio flux bin and SDSS band. We 
use the same inverted logarithmic scaling for all images. The width of the images is 24 arcsec. 
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Figure 7: Radial surface brightness profiles of the ran- 
dom super-stacks (the average of 50 random stacks per 
imaging filter), created for the coordinates of the 1 - 
2 mJy subsample using lrj masking threshold. The 
letters of the imaging filters arc indicated next to the 
curves. Error bars show the variance of the values in 
each annulus. 



to confirm the robustness of our masking and stack- 
ing method: a) jackknife analysis of the entire stacking 
process, b) analysis of the dependence of the magni- 
tude of the central peaks on the masking threshold, c) 
analysis of the dependence of the pixel value histogram 
of the stacked images on the masking threshold. 

4.1 Jackknife analysis 

To estimate the error introduced by possible outliers, 
we perform a jackknife analysis of the stacking method 
by repeating the entire processing numerous times, but 
leaving out a single cut-out image in every iteration. 

Fig. [5] shows the relative frequency of the estimated 
magnitudes of the central peaks for 709 jackknife runs 
for the g, r and i band images of the 1-2 mJy radio 
flux bin. The distributions are slightly biased towards 
the fainter magnitudes as a result of a few significant 
outliers. The few objects responsible for the faint-end 
wing of the distribution affect the stacked flux by ap- 
prox 1 per cent, thus they have a flux roughly a factor 
8 larger than the average considering a sample of 709 
objects. 

The statistical errors of the magnitudes based on the 
jackknife analysis turned out to be in the 0.1-0.2 mag 
range. In Section I5T21 we will use these error estimates 
to characterize the uncertainties of the photometry. 

4.2 Effect of the masking threshold 

Masking is required to eliminate the bright pixels of 
the cut-outs in order to avoid the contamination of 
the stacked images from moderately bright, but still 
detectable images - we are only interested in the un- 
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Figure 8: Distribution of the magnitude of the central 
peak from the jackknife stacks (histogram). The means 
of the fitted Gaussians (solid lines) are slightly offset 
to the right (see text). 



detected, faint objects. In Section 3.2.2 a masking 
threshold was introduced. The masking threshold is 
determined from the pixel value distribution of the co- 
added images by multiplying the variance of the pixel 
values with a well-chosen number. Here we show that 
setting the threshold at la is a good choice. 

Fig. IHl shows the measured magnitudes of the peaks 
detected in the centres of the stacked images for each 
SDSS band and radio luminosity bin as a function of 
the masking threshold. It is observable that the ini- 
tially increasing curves reach a plateau around la and 
then start to decrease slowly. The explanation is sim- 
ple: If the masking threshold is too low, we lose many 
pixels belonging the the investigated faint objects. On 
the other hand, if the threshold is too high the back- 
ground level of the image is increased by the pixels be- 
longing to other, nearby objects causing the contrast of 
the stacked images to be lower and the measured mag- 
nitudes to be fainter. The maxima of most of these 
curves are usually located in the 0.75er - 1.25ct range, 
so setting the masking threshold at la is obvious. 

Note that if we were to use a smoothing kernel of 
a different radius to construct the masks, the start of 
the plateau of the curves in Fig. ® would be at differ- 
ent a values. Increasing the radius of the kernel would 
move the beginning of the plateaus to the left. This be- 
haviour is a direct consequence of the masking process; 
a mask created with a larger kernel would let more pix- 
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Figure 9: Measured flux of the central peak in the 
stacked images as a function of masking threshold in 
units of a. The panels show the five SDSS bands, 
while colours refer to the radio luminosity bins (red di- 
amond: 1-2 mjy, green triangle: 2-4 mJy, blue square: 
> 4 mJy). 



els from nearby faint, point-like objects into the stack 
at a given masking threshold, because a wide kernel 
would smooth out these objects to such an extent that 
their pixel values would become less than the masking 
threshold. This is why it is important to use a smooth- 
ing kernel not much larger in diameter than the size of 
the faintest identifiable objects of the co-added images. 
This diameter is roughly the FWHM of the PSF. 

4.3 Flux deficit in the central aperture 

Because the masks are uniform in the whole cut-out 
image area, some pixels at the central source location 
can also be masked. This can introduce a slight flux 
deficit in the stack because of some marginally detected 
sources being masked out; ca. 1-5 per cent of the 
central pixels arc lost, depending on the filter (u is the 
lowest, i is the highest). If we did not mask in the 
centre, some contaminating outlier pixels would eas- 
ily get into the stack, causing a significant bias in the 
photometry (because we would overestimate the flux of 
the stacked source). So, this is a necessary trade-off in 
order to ensure the photometric quality of the stacked 
images. 

4.4 Effect of possible outlier pixels 

It is important to show that the signal emerging from 
the centres of the cut-outs after averaging comes from 
real objects and not just from a few outlier pixels. Cut- 
out images after masking will contain sky pixels (back- 
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Figure 10: Left column: Distributions of pixel values 
of the r-band cut-outs in the 1-2 mJy radio luminos- 
ity bin inside (solid curves) and outside (grey area) 
of a central circular aperture with a radius of 2 arc- 
sec for three different values of the masking threshold 
(indicated in each panel). Panels in the right column 
show the difference of the two distributions: the distri- 
butions of the background pixels subtracted from the 
distributions of the pixels within the apertures centred 
on the radio objects. 



ground), pixels of nearby faint objects that are above 
the background level, and pixels near the centres of 
the cut-outs actually belonging to the radio sources of 
our interest. The pixel values of these later ones are 
very close to the background noise level. In Fig. [TO] 
(left column), we plot the distribution of pixel values 
of the r-band cut-outs in the 1-2 mJy radio luminos- 
ity bin inside and outside of a central circular aperture 
with a radius of 2 arcsec for three different values of 
the masking threshold. The solid lines (shaded areas) 
show the normalised distributions of the pixel values 
inside (outside) the central aperture. The panels in 
the right column show the difference of the distribu- 
tions: the distributions of the background pixels sub- 
tracted from the distributions of the pixels within the 
apertures centred on the radio objects. There are two 
important things to observe in Fig. [TOj a) the distri- 
bution of background pixels have a long tail towards 
bright pixels, and b) the distribution of the pixels in- 
side the aperture is slightly shifted to the right with 
respect to the background. 

The long tail of the distribution of background pixels 
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Table 3: Signal-to- noise ratios and rms background 
noise levels of the stacked images. 

comes from the relatively bright pixels of nearby faint 
objects, like the unmasked diffuse envelopes of nearby 
galaxies. The central aperture does not contain such 
pixels because of the selection criterion we used to se- 
lect only those radio sources which are separated from 
the nearest optically detected objects by at least 3 arc- 
sec. 

The small shift between the peaks of the distribu- 
tions is due to the actual signal we are looking for. 
Optical emission from the radio objects is so faint that 
the individual pixels cannot be distinguished from the 
noise of the background, but stacking the cut-outs re- 
veals the signal. 

The lack of the bright tail of the pixel value distri- 
bution of the central pixels and the small shift between 
the distributions of background pixels and central pix- 
els confirm our claim that the detected peaks in the 
stacked images come from numerous faint objects and 
not from just a few brighter (but still faint) outlier ob- 
jects. 

5 Results 

Altogether 15 stacked images were created, one for each 
SDSS filter and for each DCS subsample. In Fig. [B] 
we have already plotted the images of the stacks. The 
PSF-like point sources are clearly visible at the centres, 
and we consider them to be the optical signal com- 
ing from the FIRST radio sources. The peaks show a 
monotonically increasing intensity towards the longer 
wavelength bands, that is the sources are red in the 
optical colours. 

5.1 Signal-to-noise ratios 

The signal-to-noise ratios and rms background noise 
levels of the stacks are summarized in Table [3] These 
signal-to-noise ratios are sufficient to calculate the pho- 
tometric properties of the central peaks, except for the 
u filter. 



5.2 Photometry 

The SDSS catalogue uses a special magnitude s ystem 
designed for noisy imaging ( Lupton et al. . 19991 ). The 
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Table 4: Average and standard deviation of the E(B- 
V) values of the DCS sample, with the correction mag 
nitrides. 



definition of the so-called 'luptitudes' is based on the 
asinh function instead of the traditionally used log- 
arithm, and it contains the scaling constant b. This 
constant was calibrated for typical SDSS frames with 
objects of average brightness. When determining the 
magnitudes of the central peaks of the stacks, it would 
be appropriate to use the same magnitude system as 
SDSS. In our case, however, the value of b had to be re- 
calibrated to account for extremely faint objects, which 
would mean the recalibration of the whole SDSS mag- 
nitude system as well. For this reason, we decided to 
use the conventional Pogson magnitudes instead. 

Because the co-added Stripe 82 images are sky- 
substracted, calibrated for atmospheric extinction and 
remapped to a common zeropoint, the transformation 
of counts into magnitudes is straightforward: 



mag A = — 2.5 log 10 counts + 30 — A^ 



(3) 



where A\ is the correction for the galactic extinction, 
discussed later in Section I5~3l 

We measure the brightness of the central peaks in a 
circular aperture with a radius of 5 arcsec. The size of 
the aperture was chosen to be slightly larger than the 
double of FWHM of the central peaks (1.5 - 2 arcsec, 
depending on the band). To characterize the photo- 
metric uncertainties, we use the standard deviations 
provided by the jackknife stacks. 



5.3 Estimating the galactic extinction 

We correct for the galactic extinction after stacking. 
We could have applied an intensity rescale to each cut- 
out image prior stacking, but that would have altered 
the noise characteristics of the images. Instead, we 
correct for galactic extinction using an average value 
of E{B — V). The average extinction is calculated 
from values of the extinction in the directions of the 
DCS sourc e coordinates. The va lues of E(B — V) are 
taken from Schleeel et al. ( 19981 ). Table |4] summarizes 
the value of the galactic extinction and the correction 
magnitudes for every imaging filter. The average value 
of E(B — V) is equal for all three subsamples to two 
decimals. Note, however, that the scatter in extinction 
values is almost as large as the average extinction it- 
self. Still, this only introduces a small error in the final 
colour index estimates, as we show below. 
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Table 5: Optical magnitudes of the central peaks iden- 
tified in the stacked images for every DCS subsamplc. 
The photometric calibration is described in Section l5~2"l 
Values are corrected for average foreground galactic ex- 
tinction. * Aperture radius: 3 arcsec. 

5.4 Optical magnitudes 

The optical magnitudes of the central peaks are sum- 
marized in Table [5] for every SDSS filter and DCS sub- 
sample. The quoted magnitudes are already corrected 
for average foreground extinction, as explained in Sec- 
tion [O] 

According to these results, the average optical mag- 
nitudes of the central peaks are about 1 mag fainter 
than the practical detection limit of the Stripe 82 co- 
added catalogue. There is a slight increase in op- 
tical brightness in all optical wavelength bands to- 
wards higher radio fluxes. The error of the photom- 
etry strongly depends on the imaging filter; the noise 
becomes really significant in the ultraviolet. 

5.5 Radial profiles 

To determine the radial profiles of the stacked objects, 
we first determine the exact centres of the central peaks 
in the images by fitting two-dimensional Gaussians to 
the pixel values. Then we average the pixels in annuli 
of dr = 0.4 arcsec centred on the fitted positions of the 
peaks. 

In Fig. [11] we plot the radial profiles of the back- 
ground corrected stacks for every filter and DCS sub- 
sample. The central peaks of the profiles are clearly 
visible, just like the faint, extended components that 
decrease with radius. 

We fit the profiles with the combination of a Gaus- 
sian and an exponential function. These fits are also 
displayed in Fig. Qj] (solid lines). The fitting function 
is formulated the following way: 



I(r) 
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exp 



r 

2^ 



B exp 
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(4) 



The fitting is done in two phases. First, we fit the 
outer part (3 < r < 8 arcsec) of the profiles with the 
exponential term to obtain the value of B and the 
scalelength r$. Then, the resulting exponential fit is 
subtracted from the original profile. Finally, the resid- 
uals are fitted with Gaussians in the 0-5 arcsec radius 
range. The fitted parameters are summarized in Ta- 
ble |H1 The central PSF-like peaks cut off around r ~ 2 
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r [arcsec] r [arcsec] r [arcsec] 

Figure 11: Surface brightness profiles of the stacked images around the central peaks. The graphs in each 
column show the profiles for a given radio flux bin for all five SDSS imaging filters. Solid lines indicate the 
model fits; dashed lines are the Gaussian components; dotted lines are the exponential components. Error bars 
represent the variance of the values in each annulus. 
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Si n t,i400 [mJy] 
1-2 2-4 > 4 

~A 9.56 ± 5.90 19.96 ± 4.92 21.62 ± 3.20 

a 0.51 ± 0.32 0.62 ± 0.16 0.63 ± 0.10 

B 1.57 ± 0.82 2.55 ± 1.59 3.97 ± 2.23 

rg 28.70 ± 56.13 7.24 ± 5.41 9.38 ± 9.09 

A 44.23 ± 2.47 39.59 ± 2.79 46.26 ± 2.82 

a 0.87 ± 0.05 0.68 ± 0.05 0.75 ± 0.05 

B 1.59 ± 0.90 4.47 ± 1.46 4.61 ± 1.26 

rp 6.01 ± 4.14 4.63 ± 1.36 4.54 ± 1.11 

A 87.24 ± 2.78 63.10 ± 3.84 87.39 ± 4.56 

a 0.82 ± 0.03 0.62 ± 0.04 0.72 ± 0.04 

B 2.82 ± 1.04 10.00 ± 3.41 9.16 ± 2.05 

r 5.98 ± 2.82 3.52 ± 0.89 3.85 ± 0.67 

A 173.73 ± 6.60 139.96 ± 7.76 147.14 ± 8.10 

a 0.82 ± 0.03 0.72 ± 0.04 0.67 ± 0.04 

B 4.39 ± 2.13 9.71 ± 3.04 20.09 ± 3.47 

rp 6.03 ± 3.94 5.34 ± 1.68 3.55 ± 0.44 

A 394.57 ± 26.77 279.43 ± 27.24 405.69 ± 20.20 

cr 0.74 ± 0.05 0.62 ± 0.06 0.71 ± 0.04 

B 4.27 ± 2.13 18.69 ± 8.40 25.27 ± 11.34 

r 18.12 ± 36.27 10.37 ± 8.10 3.48 ± 1.31 



Filter 


<Sint,i400 [mJy 






1-2 2-4 


> 4 


g-r 


0.7 ± 0.2 0.6 ± 0.2 


0.7 ± 0.1 


r — i 


0.6 ± 0.2 0.6 ± 0.2 


0.6 ± 0.1 


i — z 


0.8 ± 0.2 0.6 ± 0.2 


0.6 ± 0.1 



Table 7: Optical colour indices of the central peaks 
identified in the stacked images for every DCS sub- 
sample. 



5int,1400 




bin [mJy] 




1 - 2 


-2.9 ±0.3 


2-4 


-2.5 ±0.1 


>4 


-2.2 ±0.3 



Table 8: Average spectral indices of the central peaks. 



Table 6: Fitted parameters of the radial profiles: A 
[nJy arcsec -1 ], a [arcsec], B [nJy arcsec -2 ], ro [arcsec]. 



- 3 arcsec, but the exponential component is traceable 
to ca. r ~ 15 arcsec. The u profiles in the 1-2 mJy and 
> 4 mJy bins do not have good enough data quality to 
really trust the fitted values. 

5.6 Possible origins of the exponential 
component 

The existence of an exponentially decreasing compo- 
nent of the radial profiles is obvious from Fig. [TTJ The 
value of the scale parameter ro are listed in Table [5) 
Excluding the low quality fits of the u filter, the value 
of ro is consistently around r ~ 3.5 - 5.5 arcsec. 

It is very important to consider that the PSF of 
SDSS images already shows an extended component. 
To investigate the contribution of these extended wings 
to the radial profiles of our stacks, we applied our stack- 
ing algorithm to images of faint stars. The stacks of 
stars showed a similar exponential tail with a scale pa- 
rameter ro in the range of 2 - 3 arcsec, depending on 
the photometric band. These values are significantly 
smaller than the scale lengths derived from the stacks 
of our radio-selected objects. The radial profiles of 
our stacks clearly show an excess over the stellar PSF 
at the tails. For more discussion on the effect of the 
PSF wings on the r adial profiles of stack ed objects, see 
Zibetti et all (|2004[) and lde Jond <|2008ft . 

It is very risky to draw any conclusions from the ob- 
servation of the exponential component, especially in 
the absence of information about the redshift of the 
radio-selected objects. Because the redshifts are un- 
known, we cannot rescale our sources to the same phys- 
ical scale before the stacking. The physical scales could 
therefore be very different between the objects that are 
stacked together. 

However, it is interesting to see what the measured 



angular scalelength translates to at different redshifts. 
If we make the assumption that the majority of the 
DCS sample lies around z ~ 0.1, the corresponding 
scalelength is 6-10 kpc. At z ~ 1, the scalelength would 
be 28-45 kpc. The former range is consistent with the 
size of individual galaxies, which might suggest that 
we see the host galaxies of the quasars. The latter 
scalelength could correspond to the diameter of merg- 
ing systems. It can also be imagined that the excess in 
surface brightness at large radii comes from undetected 
galaxies clustered around the radio sources. 

5.7 Spectral energy distributions 

One important finding is that the stacked objects have 
unusually red average colours, as listed in Table [7] The 
colour indices seem to be independent from the appar- 
ent radio luminosity. In Fig. [12j we plot the SEDs of 
the central peaks of the stacked images for all three 
DCS subsamples (solid lines with symbols). The op- 
tical spectral indices {a v ) are determined by fitting 
power-law functions to the fluxes, and the results are 
listed in Table ® 

Because the redshift distribution of the DCS sample 
is unknown, the colours of the stacked object might 
not correctly represent the average spectral slope of 
the sample. In order to obtain the correct colours, one 
should perform k-correction for each DCS source indi- 
vidually - obviously, this cannot be done in this study. 
However, if we assume that the majority of the sources 
have a power-law SED in the optical and ultraviolet 
regime, with similar spectral slopes, the k-correction 
magnitudes in this case are independent of wavelength, 
and hence the colour indices remain unchanged. 

However, we note that if the objects had a redshift 
distribution similar t o that of the LBD S sam ple men- 
tioned in Section [T31(jWaddington et all 200lh . that is, 
the majority of the objects had z < 1.5, based on the 
average spectral index of a v = —2.5, then the average 
k-correction would be ~ 0.9 mag in every band. 
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Figure 12: Optical SEDs (symbols connected with solid 
lines) of the central peaks identified in the stacked im- 
ages for every DCS subsample. The short grey lines 
indicate the fluxes of the individual objects of the cross- 
check sample described in Section H3 The dashed line 
is the average SED of the cross-check sample. 



Vanden Berk et al.l ([20011 ) calculated composite 



spectra of more than 2200 SDSS quasars and they 
found a v = —0.44 for the spectral continuum. 



Ilvezic et al. ( 2002) determined spectral indices individ- 
ually for 6868 quasar spectra with a mean of a v — 
—0.45. The bulk of th eir distributi o n is in the range 

" (|2002h suggested 



Gregg et al 



of -1 < a v < 0.2 
ol v < — 1 as a definition of a red quasar. They found 
two sources with a v ~ —3.7 a nd ay < —4.6, which 



they called extraordinarily red. iRichards et al.l ([2006) 
defined Type 1 quasars by a — — 0.5±0.3 and reddened 
Type 1 quasars by ol v < —1, which both have broad 
emission lines. The values for our stacked sources are 
in the range —2.9 < a v < —2.2. However, note that 
we cannot account for the effect of emission lines when 
calculating the spectral indices, and therefore we have 
a v values from fitting of the entire spectrum (contin- 
uum and lines). Measurements of a by others usually 
based on the spectral continuum only. 

If we account the red colours to intrinsic extinction, 
the values of— 2.9 < a v < —2.2 require unusually high 
column densities of the obscuring dust. The absolute 
value of the spectral index becomes slightly smaller 
with increasing radio flux. This might suggest - if we 
assume that the redshift distribution of the objects is 
narrow, or at least similar for all radio flux bins - that 
the brighter radio sources are slightly less reddened on 
average. 



6 Discussion 

Because the individual objects of the DCS sample are 
undetectable in the optical images, we can only discuss 
the average optical properties of the objects based on 
their stacked images. We hope that these average val- 



ues represent the original distribution of the sources 
well. With the lack of direct detections, the homogene- 
ity of the sample cannot be proven. However, it is still 
worth comparing the DCS sample with samples of very 
similar selection criteria, except that the radio-selected 
objects are also detected in the optical. 

6.1 Construction of the cross-check 
sample 

In order to compare the average optical properties of 
the DCS objects to better-known objects, we construct 
a cross-check sample of radio selected sources, but in 
this case with identified optical counterparts. We select 
FIRST radio point sources that have optical counter- 
parts within 1.5 arcsec in the Stripe 82 catalogue, but 
not in the singly observed (thus less faint) SDSS Data 
Release 6 (DR6) catalogue within 3 arcsec. The 1.5- 
arcsec matc h ing ra dius was chosen in accordance with 
Ivezic et alJ (|2002l) . who matched FIRST with SDSS 
using the same separation limit. Based on their study, 
we expect similar completeness (~ 85 per cent), and a 
slightly higher level of contamination than their 3 per 
cent, because the larger source density in the Stripe 82 
co-added catalogue increases the false-positive match 
rate. All other selection criteria for the FIRST sources 
are the same as described in Section [5J 

Because the co-added Stripe 82 is a significantly 
deeper survey than the singly observed SDSS DR6, 
we expect to find objects that are too faint to be de- 
tected in DR6 but are visible in the co-added Stripe 
82 images. Indeed, 1349 FIRST sources were found 
with optical magnitudes in the range 24 > > 22 
mag. We apply further restrictions on the photomet- 
ric quality of the cross-check sample (i.e. the error in 
u and g model magnitudes must be less than 0.8 and 
0.5 mag, respectively). The final sample consists of 394 
co-added Stripe 82 sources, which are optical counter- 
parts of isolated, compact FIRST radio sources with 
integrated luminosities larger than iSi n t,i400 > 1 m Jy- 

We have published a catalogue on the multiwave- 
length photometric properties of the cross-check ob- 
jects, compiled from the Stripe 82 co-added sur- 
vey (optical), the FIRST survey (radio) and two in- 
frared surveys: the United Kingdom Infrared Tele- 
scope (UKIRT) Infrared Deep Sky Survey (UKIDSS) 
and the Wide-field Infrared Survey Explorer (WISE). 
The infrared cross-identification is described in Sec- 
tion 16.31 The catalogue is available on-line at 



http : //www . vo . elte . hu/ doublestacking 

We plot the optical fluxes of the cross-check sample 
in Fig. rT2] as grey dashes. The average fluxes of the 
cross-check sample (crosses connected with a dashed 
line) and the fluxes of the stacked objects (points con- 
nected with solid lines) are also plotted in the same 
plot. The average magnitudes calculated from the 
stacked images are about 1.2—1.6 mag fainter than 
those of the cross-check sample, but the slopes of the 
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Figure 13: Colour-colour diagrams of the stacked objects (filled circles - red: 1-2 mJy, green: 2-4 mJy, blue: 
> 4 mJy), the cross-check sample (orange crosses), normal Type 1 quasars (black dots), red Type Is (purple 
squares), and Type 2s (red contour lines). See the description of the samples in Section 16721 A colour version 
of this plot is available on-line. 




Figure 14: Optical and infrared SEDs of the DCS sample (coloured dots: red: 1-2 mJy, green: 2-4 mJy, 
blue: > 4 mJy subsamples) and the cross-check sample (black crosses). The data points in all infrared bands 
are individual detections, while the optical fluxes of the stacked objects are average values. We plot only those 
objects of the cross-check sample that have matching detections in the K band. The average optical fluxes of 
the entire cross-check sample are plotted with black squares. (A colour version of this figure is available on-line.) 
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SEDs are similar. 

6.2 Comparison with the cross-check 
sample and optically detected 
quasars 

To evaluate the average optical properties of the DCS 
objects we compare them to three samples of quasars. 
The first sample is the fifth edition of the SPSS quasar 



catal ogue containing 105783 objects (jSchneider et al 

The catalogue is composed of SDSS quasars 
with luminosities higher than Mj < —22.0, having 
at least one emission line with FWHM higher than 
1000 km s _1 or showing interesting/complex absorp- 
tion features. The quasar redshifts range from 0.065 
to 5.46, with a median value of 1.49. The major- 
ity of the sample can be classified as normal (unred- 
dened) Type 1 quasars. The s econd quasar sampl e 
consists of 51 red Type 1 quasars ( Urrutia et all l2009h . 
lUrrutia et aH matched the F IRST-2MASS red quasar 
survey (jGlikman et~all . l2007l) with SDSS, and compiled 
this sample of spectroscopically confirmed red Type 1 
quasars. The third sample contains 887 spectroscopi- 
cally confirmed Type 2 quasars fro m SDSS with red- 



shifts z < 0.83 (jReves et all [2008). The stacked ob 



jects are about 1.4 mag fainter in the i band than the 
average of the cross-check sample and the latter is 3 
mag fainter than the red Type 1 sample. 

In Fig.[T3]we plot the optical colour-colour diagrams 
of the stacked objects, the cross-check sample and the 
three quasar samples. It is clear that on average our 
DCS sample is significantly redder than normal Type 1 
quasars. The r — i and i — z colours suggest that 
the DCS sample is more likely to consist of reddened 
Type Is than Type 2s. 

6.3 Infrared counterparts 

If we accept the assumption that our radio sources 
are AGNs heavily obscured by dust surrounding the 
central regions, it is straightforward to look for ex- 
cess in their infrared emission. For this reason, we 
looked for infrared counterparts to the DCS sample 
in the extragalactic near-infra red UKIDSS Lar ge Area 
Survey (LAS) DR7 catalogue (|Dve et all [20061) and in 
the WI SE preliminary relea se mid-infrared source cat- 



alogue (jWright et all l2010h . Within the sky regions 



overlapping with SDSS Stripe 82, about 4 per cent of 
the radio sources (72 objects) were detected by the 
UKIDDS in any of the infrared bands, using a match- 
ing radius of 1.5 arcsec. Only a few dozen objects have 
matching counterparts in the WISE catalogue in at 
least one band. This detection rate is about 5 per cent. 
The infrared magnitudes are close to the detection lim- 
its of both surveys and the scatter of the magnitudes 
is small. This means that only the very bright end of 
the sample is detected in infrared. 

The cross-check sample has a much higher detection 



rate, with 123 sources detected in the near-infrared (37 
per cent) and 33 in the mid-infrared (32 per cent). The 
distribution of the infrared magnitudes is similar to 
that of the DCS sample, with cut-off at the faint end 
because of the magnitude limit. In Fig.[T4l we plot the 
optical to mid-infrared SED of our samples: the indi- 
vidually detected objects of the DCS sample (colour 
dots) and the cross-check sample (black crosses). The 
data points in all infrared bands are individual detec- 
tions, while the optical fluxes of the stacked objects 
are average values for the three DCS subsamples. We 
plot the average optical fluxes of the cross-check sample 
with black squares. 

The cross-check sample has an average SED with a 
slope of a v ~ —2.0 (see Fig. [TJ]), less steep than the 
slope of the average optical SED of the DCS sample. 
Exact physical interpretation is difficult, however, as 
we do not have information on the composition of the 
samples. If we assume that the objects in the samples 
are intrinsically similar, the gradual steepening of the 
spectral slopes with decreasing optical brightness can 
be attributed to obscuration by dust. 



Glikman et al.1 (|2004|) constructed a sample of bright 



near-infrared sources that are detected at radio wave- 
lengths but unidentified on the Palomar Observatory 
Sky Survey (POSS) plates, with the aim of finding 
dust-obscured quasars. They have defined a region in 
the R—K, J — K colour plane, in which 50 per cent of 
the radio-selected objects are highly reddened quasars 
(J — K > 1.7 and R — K > 4). To compare our sam- 
ples to these criteria, we converted SPSS magnitudes t o 
Johnson R using the equations of I Jester et al. (2005). 
There are 17 sources of our cross-check sample that 
have detections in both J and K infrared bands. For 
objects not detected in the J band (the majority), we 
assume a lower limit of the J — K colour index using 
the limiting magnitude of the survey. In Fig. 1151 we 
plot the cross-check objects in the R — K, J — K colour 
plane. Only a few objects have valid J magnitudes (di- 
amonds); object with lower limits in J are plotted with 
grey arrows. The majority of our sources (77 per cent) 
are in the red quasar region of the plane. 

Only a few objects of the PCS sample were detected 
in infrared individually. It is very likely that the op- 
tical flux distribution of these objects differs signifi- 
cantly from the distribution of the entire PCS sample. 
Because the PCS sources are undetected in the opti- 
cal, we can easily determine an upper limit on their R 
fluxes; they should be fainter than the Stripe 82 co- 
added survey limit (R > 24.2). 

Therefore, we are able to estimate a lower limit on 
the R — K colour of the PCS sources with infrared 
detections. We also calculate a lower limit on J — K 
in the absence of reliable J band fluxes. The results 
are plotted in Fig. [15] with coloured symbols (arrows) , 
one for each PCS subsample. It is apparent that the 
R — K values are in the red quasar region, even if we 
assume brighter R fluxes. The J — K colours are also 
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Figure 15: Colour indices of the cross-check sample 
and the DCS sample in the R — K; J — K plane. Di- 
amonds denote cross-check sample objects with valid 
J band data, arrows denote cross-check objects with 
the J band lower limit only and coloured arrows de- 
note stacked objects (red: 1-2 mjy, green: 2-4 
mjy, blue: > 4 mjy subsamples). The arrows mark 
the lower limits on R — K and J — K . The dashed lines 
bound the upper right region, where 50 per cent of the 
source s are red quasars, according to Idikman et al 
(2003). 



in that region, but much closer to the border. 

It is hard to estimate the possible fraction of red 
quasars in our samples based on the available infrared 
data. As for the cross-check sample, with a 34 per 
cent detection rate and with at least 77 per cent of 
the infrared-matched objects lying in the red quasar 
region of the R — K, J — K plane, we find that at 
least 13 per cent of all objects are undoubtedly dust- 
obscured quasars. However, this is a weak lower limit 
on their number, because we only see the brightest in- 
frared counterparts as a result of the relatively high 
infrared flux limit of the observations. 



7 Summary 

We have presented an analysis based on an image stack- 
ing technique, to reveal the visible wavelength light 
from the unresolved sources of the FIRST radio survey 
that remain undetected in the SDSS Stripe 82 co-added 
catalogue. The sample consists of 2116 objects, which 
have been divided into three subsamples according to 
the radio flux. We stacked cut-out images centred on 
the object coordinates from the SDSS Stripe 82 co- 
added survey. The main results of our study are as 
follows. 

1. We have presented an image stacking algorithm 
with an efficient masking technique to reveal the 
average optical emission of 2116 optically unde- 
tected sources (DCS sample) selected based on the 



VLA FIRST radio survey. 

2. We have detected the average optical emission of 
these radio sources (with 26.6 mag 5a detection 
limit in r band), and we have calculated the ugriz 
magnitudes of the peaks apparent in the centres 
of the stacked images (see Table [5]). The colours 
of these peaks imply a steep, red optical spectrum 
with spectral indices of —2.9 < a v < —2.2. 

3. The average surface brightness profiles of the 
stacked objects show Gaussian peaks in the cen- 
tres, which continue in outwards fading exponen- 
tial components traceable to ~ 15 arcsec. The 
surface brightness of these components peaks at 
approximately 2-3 mag arcsec -2 fainter than the 
peak surface brightness of the central point-like 
components. 

4. We have created a radio-selected sample with faint 
optical detections in the SDSS Stripe 82 catalogue 
{cross-check sample), similar to the DCS sam- 
ple. We have identified infrared counterparts to 
the cross-check objects from the UKIDSS near- 
infrared and WISE mid-infrared source catalogue 
(with 37 per cent detection rate), and we have 
composed a catalogue of the optical, infrared and 
radio properties of the sample. 

5. We have compared the optical colour indices of the 
radio-selected and cross-check objects with var- 
ious spectroscopically identified quasar samples. 
We have concluded that the distribution of the 
colours are similar to that of dust-reddened Type 1 
quasars, although there is large scatter in the data. 

6. We have identified the infrared counterparts of 
the DCS sources with a 5 per cent detection rate. 
We have investigated the distribution of the cross- 
check and DCS sources in the R— K, J — K colour 
plane, and we have found that the majority of the 
objects lie in a region where ~ 50 per cent of the 
sources are red quasars. 

7. Consequently, the sources of the DCS sample and 
the cross-check sample are very good red quasar 
candidates, suitable for future deeper observations 
and quasar studies. We assume that the forthcom- 
ing optical and infrared survey telescopes, such as 
the Large Synoptic Survey Telescope (LSST) and 
the James Webb Space Telescope (JWST), might 
discover many more obscured quasars than previ- 
ously expected. 
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